1 Setup R Environment

library(knitr)
library(rmdformats)
## Global options
options(max.print="75")
opts_chunk$set(#echo=FALSE,
               cache=TRUE,
               prompt=FALSE,
               tidy=TRUE,
               comment=NA,
               message=FALSE,
               warning=FALSE, 
               results='asis'
               )
opts_knit$set(width=75)
print("notebook options imported successfully!")
## [1] "notebook options imported successfully!"
## Load Required Libraries & Set Paths
library(tidyverse)
library(tidyquant)
library(lubridate)
library(timetk)
library(DT)
library(PerformanceAnalytics)
library(ROI.plugin.quadprog)
library(ROI.plugin.glpk)
# library(ROI.plugin.symphony)
setwd("~/Documents/GitHub/rob_gordon_2018")

2 Get Stock Data

stocks <- c("SPY", "IYR", "QQQ", "TLT", "GLD", "IWM", "EEM", "DBC", "EFA")
stock_data <- tq_get(stocks
                     , get  = "stock.prices"
                     , from = Sys.Date() - months(12)
                     , to = Sys.Date()) 
stock_data %>% head(10)
symbol date open high low close volume adjusted
SPY 2017-11-15 256.62 257.22 255.63 256.44 80811500 251.8421
SPY 2017-11-16 257.52 259.04 257.47 258.62 67777000 253.9830
SPY 2017-11-17 258.22 258.59 257.77 257.86 75756800 253.2366
SPY 2017-11-20 258.14 258.52 257.86 258.30 48075500 253.6688
SPY 2017-11-21 259.18 260.20 258.26 259.99 69176800 255.3284
SPY 2017-11-22 260.00 260.15 259.57 259.76 45033400 255.1026
SPY 2017-11-24 260.32 260.48 260.16 260.36 27856500 255.6918
SPY 2017-11-27 260.41 260.75 260.00 260.23 52274900 255.5642
SPY 2017-11-28 260.76 262.90 260.65 262.87 98971700 258.1568
SPY 2017-11-29 263.02 263.63 262.20 262.71 77512100 257.9997

3 Visualize Stock Prices

  • Kept only 4 symbols in order to fit
  • Will include all in dashboard via interactive filtering
# using ggplot2
stock_data %>%
    filter(symbol %in% c("SPY", "GLD", "QQQ", "DBC")) %>%
    ggplot(aes(x = date, y = adjusted, color = symbol)) +
    geom_line(size = 1) +
    geom_bbands(aes(high = high, low = low, close = close), ma_fun = SMA, n = 5,show.legend=TRUE) + 
    labs(title = "Daily Stock Prices",
         x = "", y = "Adjusted Prices", color = "") +
    facet_wrap(~ symbol, ncol = 2, scales = "free") +
    scale_y_continuous(labels = scales::dollar) +
    theme_tq() + 
    scale_color_tq()

4 Calculate Returns

4.1 Simple Returns Data

mo_returns <- stock_data %>% 
  group_by(symbol) %>%
  tq_transmute(select = adjusted, 
                 mutate_fun = periodReturn, 
                 period     = "monthly", 
                 col_rename = "returns") 

mo_returns %>% 
  mutate(returns = round(returns, 6)) %>% 
  arrange(desc(symbol), desc(date)) %>% datatable(class="compact", options = list(scrollX = TRUE, dom = "tp"))
# simple returns = Rt = [ Pt - Pt-1 ]/ Pt-1
# simple_returns <- 
#   stock_data %>% 
#   group_by(symbol, year = year(date)) %>%
#   summarise(adjusted = last(adjusted, order_by = year)) %>% 
#   mutate(return = (adjusted - lag(adjusted, order_by = year)) / lag(adjusted, order_by = year)
#          , return = round(return, 6)) %>% 
#   arrange(desc(symbol), desc(year)) 
#   # %>% filter(symbol == "SPY")
# datatable(simple_returns, class="compact", options = list(scrollX = TRUE, dom = "tp"))

4.2 Chart: Returns

  • stocks filtered to fit
mo_returns %>% 
  filter(symbol %in% c("SPY", "GLD", "QQQ", "DBC")) %>%
    ggplot(aes(x = date, y = returns, fill = symbol)) +
    geom_line(aes(color=symbol)) +
    geom_hline(yintercept = 0, color = palette_light()[[1]]) +
    scale_y_continuous(labels = scales::percent) +
    labs(title = " Monthly Returns",
         subtitle = "stocks limited to fit",
         y = "Monthly Returns", x = "") + 
    facet_wrap(~ symbol, ncol = 2,scales = "free") +
    theme_tq() + 
    scale_fill_tq()

4.3 Chart: Smoothed Returns

mo_returns %>% 
  filter(symbol %in% c("SPY", "GLD", "QQQ", "DBC")) %>%
    ggplot(aes(x = date, y = returns, fill = symbol)) +
    geom_smooth(aes(color=symbol),method = 'loess' , formula = y ~ x, se=FALSE) +
    geom_hline(yintercept = 0, color = palette_light()[[1]]) +
    scale_y_continuous(labels = scales::percent) +
    labs(title = " Monthly Returns",
         subtitle = "stocks limited to fit",
         y = "Monthly Returns", x = "") + 
    facet_wrap(~ symbol, ncol = 2,scales = "free") +
    theme_tq() + 
    scale_fill_tq()

5 Growth of $1000 (6mo)

5.1 Data Table

init.investment <- 1000
growth <- mo_returns %>% arrange(date) %>%
  mutate(final_value = init.investment * cumprod(1 + returns)) %>%
  arrange(desc(final_value)) 
growth %>% filter(date == max(date)) %>% select(-date)
symbol returns final_value
QQQ -0.0272053 1090.9494
SPY -0.0015889 1072.8945
IWM -0.0038659 1038.9866
IYR 0.0266291 1016.4514
DBC -0.0536557 996.2755
GLD -0.0044290 944.2385
EFA 0.0036824 941.7529
TLT 0.0069282 925.4958
EEM 0.0196629 893.9512

5.2 Visualization

growth %>% ggplot(aes(x = date, y = final_value, color = symbol)) +
    geom_line() +
    # geom_smooth(method = "loess") +
    labs(title = "Individual Portfolio: Comparing the Growth of $1K",
         subtitle = "Quickly visualize performance",
         x = "", y = "Investment Value") +
    theme_tq() + theme(legend.position = "right") +
    scale_y_continuous(labels = scales::dollar)

## winners (top 5)

winners <- growth %>% ungroup() %>% filter(date == max(date)) %>% mutate(rank = row_number()) %>% 
    top_n(5, final_value) %>% select(rank, symbol, final_value)

winners %>% arrange(desc(final_value)) %>% ggplot(aes(x = symbol, y = final_value, 
    fill = symbol)) + geom_bar(stat = "identity") + # geom_smooth(method = 'loess') +
labs(title = "Individual Portfolio: Comparing the Growth of $1K", subtitle = "Quickly visualize performance", 
    x = "", y = "Investment Value") + theme_tq() + theme(legend.position = "right") + 
    scale_y_continuous(labels = scales::dollar) + geom_text(aes(label = round(final_value, 
    2), vjust = 3, size = 10))

winners$symbol

[1] “QQQ” “SPY” “IWM” “IYR” “DBC”

6 Portfolio Optimization

6.1 Create Time-Series Object for Optimization

# convert to ts object
stock_returns <- mo_returns %>% filter(symbol %in% winners$symbol) %>% spread(key = symbol, 
    value = returns) %>% tk_xts(date_var = date) %>% na.omit()

# stock_returns %>% class()

7 Optimize & BACKTEST

7.1 Get Historical Stock Returns (3 Years)

hist_returns <- tq_get(winners$symbol, get = "stock.prices", from = Sys.Date() - 
    months(36), to = Sys.Date()) %>% group_by(symbol) %>% tq_transmute(select = adjusted, 
    mutate_fun = periodReturn, period = "daily", type = "arithmetic", col_rename = "returns")

hist_returns_ts <- hist_returns %>% spread(key = symbol, value = returns) %>% 
    tk_xts(date_var = date) %>% na.omit()

7.2 Run Minimum Variance Optimization

  • set objective: risk minimization (variance)
  • contraint 1: invest all funds
  • contrstaint 2: long only
# Markowitz Optimization involves minimizing the risk(variance) of returns and maximizing returns.
# https://rdrr.io/cran/PortfolioAnalytics/f/inst/doc/ROI_vignette.pdf
library(PortfolioAnalytics)

# Loading all suggested packages to fix ROI issue
# list.of.packages <- c("foreach", "DEoptim", "iterators", "fGarch", "Rglpk", "quadprog", "ROI", "ROI.plugin.glpk", "ROI.plugin.quadprog", "ROI.plugin.symphony", "pso", "GenSA", "corpcor", "testthat", "nloptr", "MASS", "robustbase")
# new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
# if(length(new.packages)) install.packages(new.packages)

#' Set up initial portfolio with basic constraints.
init.portf <- portfolio.spec(assets = colnames(hist_returns_ts))
# Add full investment constraint to the portfolio object
init.portf <- add.constraint(portfolio=init.portf, type="full_investment")
# Add long only constraints
init.portf <- add.constraint(portfolio=init.portf, type="box", min_sum=0.99, max_sum=1.01)

# Add objective to minimize variance
portf_minvar <- add.objective(portfolio=init.portf, type="risk", name="var")
# Dummy objectives for plotting and/or further analysis
portf_minvar <- add.objective(portf_minvar, type="return", name="mean", multiplier=0)
portf_minvar <- add.objective(portf_minvar, type="risk", name="ES", multiplier=0)
portf_minvar <- add.objective(portf_minvar, type="risk", name="StdDev", multiplier=0)

opt_gmv <- optimize.portfolio(R=hist_returns_ts,
                              portfolio=portf_minvar,
                              optimize_method="ROI",
                              trace=TRUE,
                              search_size=5000)

### Print Min Variance Achieved via Optimization ###
extractStats(opt_gmv)[1]
    ES

0.01513927

### Print Weights ###
weights <- extractWeights(opt_gmv)
weights_tbl <- as.tibble(as.list(weights)) %>% gather("symbol", "weight", 1:length(weights)) %>% arrange(desc(weight))
weights_tbl
symbol weight
IYR 0.4202591
DBC 0.3657620
SPY 0.2139789
IWM 0.0000000
QQQ 0.0000000 
### Bar Chart Weights ###
opt_gmv %>% chart.Weights(plot.type="barplot", ylim=c(0,1), las=1, colorset = rainbow12equal, main="Optimal Weights", cex.axis = .8 )

### Pie Chart Weights ###
library(scales)
pie <- ggplot(weights_tbl, aes(x="", y=weight, fill=symbol)) +
  geom_bar(width = 1, stat = "identity") + 
  coord_polar("y", start=0) +
  scale_fill_brewer("Blues") + #blank_theme() +
  theme(axis.text.x=element_blank())+
  geom_text(aes(y = weight/3 + c(0, cumsum(weight)[-length(weight)]), 
                label = percent(weight)), size=5)
pie

paste("Min Variance Achieved: ", round(opt_gmv$objective_measures$ES,5))

[1] “Min Variance Achieved: 0.01514”

### Chart Risk Reward ### 
chart.RiskReward(opt_gmv, rp=TRUE, risk.col = "ES", chart.assets=TRUE)

annualized.moments <- function(R, scale=12, portfolio=NULL){
    out <- list()
    out$mu <-    matrix(colMeans(R), ncol=1)
    out$sigma <- cov(R)/scale
    return(out)
  }

### CALCULATE & PLOT EFFICIENT FRONTIER ### 
  prt_ef <- create.EfficientFrontier(R=12*hist_returns_ts, portfolio=portf_minvar, type="mean-StdDev", 
                                      match.col = "StdDev", momentFUN="annualized.moments", scale=12)
  ### Chart Efficient Frontier ### 
  xlim <- range(prt_ef$frontier[,2])*c(1, 1.5)
  ylim <- range(prt_ef$frontier[,1])*c(.80, 1.05)
  chart.EfficientFrontier(prt_ef, match.col="StdDev", chart.assets = TRUE, n.portfolios = 1000
                          , labels.assets = TRUE, xlim=xlim, ylim=ylim, main = "Efficient Frontier")
  points(with(annualized.moments(12*stock_returns, scale=12), cbind(sqrt(diag(sigma)), mu)), pch=19 ) 
  text(with(annualized.moments(12*stock_returns, scale=12), cbind(sqrt(diag(sigma)), mu)), 
       labels=colnames(stock_returns), cex=.8, pos=4) 

  ### Chart EF Weights ###
  chart.EF.Weights(prt_ef, match.col="StdDev", colorset = rainbow12equal, main = "Efficient Frontier")

8 Visualize Portfolio Optimization

8.1 Perform Backtest

# https://rdrr.io/cran/PortfolioAnalytics/f/inst/doc/ROI_vignette.pdf

bt_gmv <- optimize.portfolio.rebalancing(R=hist_returns_ts, portfolio=portf_minvar,
                                         optimize_method="ROI",
                                         trace = TRUE,
                                         rebalance_on="months",
                                         training_period=36,
                                         # trailing_period=trailing, 
                                         rolling_window = 12
                                         )

summary(bt_gmv)$annualized_returns[[1]]

[1] 0.1133622

summary(bt_gmv)$annualized_StdDev[[1]]

[1] 0.1018277

bt_weights <- summary(bt_gmv)$weights
bt_ret <- summary(bt_gmv)$portfolio_returns

charts.PerformanceSummary(cbind(bt_ret, opt_gmv$R), main="Optimization Performance", method = "StdDev", wealth.index = TRUE, begin = "axis")

x <- Return.portfolio(hist_returns_ts, bt_weights, wealth.index = TRUE, verbose = TRUE)
chart.CumReturns(x$returns, col=bluemono, begin="axis")

chart.StackedBar(x$BOP.Weight)

chart.StackedBar(x$BOP.Value)

# CHART BACKTEST WEIGHTS
bt_gmv %>% chart.Weights(ylim=c(0,1), las=1, main="Optimal Weights", cex.axis = 3 ,cex.legend = 3, col=rainbow12equal) # , legend.loc = "topright"

# ret.bt.gmv <- do.call(cbind, lapply(bt_gmv, function(x) summary(bt_gmv)$portfolio_returns))
# colnames(ret.bt.gmv) <- c("min ES", "min ES RB", "min ES Eq RB", "stuff", "other stuff")

Rb <- Return.portfolio(hist_returns_ts, weights = c(.2,.2,.2,.2,.2), wealth.index = TRUE, verbose = TRUE)
chart.RelativePerformance(as.xts(x$returns), as.xts(Rb$returns), colorset = rainbow12equal, main = "Performance Relative to Even Dist")

9 Get Performance

9.1 descriptive statistics

symbol ArithmeticMean GeometricMean Kurtosis LCLMean(0.95) Maximum Median Minimum NAs Observations Quartile1 Quartile3 SEMean Skewness Stdev UCLMean(0.95) Variance
SPY 0.0060 0.0054 0.0073 -0.0146 0.0564 0.0059 -0.0691 0 13 -0.0016 0.0319 0.0094 -0.7344 0.0340 0.0265 0.0012
IYR 0.0017 0.0013 -0.4254 -0.0176 0.0406 0.0023 -0.0666 0 13 -0.0239 0.0266 0.0089 -0.6159 0.0319 0.0210 0.0010
QQQ 0.0077 0.0067 -0.0764 -0.0197 0.0876 0.0060 -0.0860 0 13 -0.0129 0.0280 0.0126 -0.2076 0.0454 0.0351 0.0021
TLT -0.0057 -0.0059 -1.5196 -0.0190 0.0286 -0.0114 -0.0326 0 13 -0.0286 0.0131 0.0061 0.1430 0.0220 0.0076 0.0005
GLD -0.0042 -0.0044 -0.7278 -0.0162 0.0323 -0.0066 -0.0361 0 13 -0.0208 0.0063 0.0055 0.3733 0.0199 0.0078 0.0004
IWM 0.0039 0.0029 1.3954 -0.0230 0.0616 0.0098 -0.1099 0 13 -0.0040 0.0256 0.0124 -1.1468 0.0446 0.0308 0.0020
EEM -0.0076 -0.0086 -0.4608 -0.0354 0.0830 -0.0058 -0.0876 0 13 -0.0377 0.0197 0.0128 0.1761 0.0461 0.0202 0.0021
DBC 0.0002 -0.0003 -1.1508 -0.0198 0.0342 0.0075 -0.0562 0 13 -0.0243 0.0278 0.0092 -0.5503 0.0331 0.0202 0.0011
EFA -0.0040 -0.0046 0.2145 -0.0249 0.0502 0.0037 -0.0813 0 13 -0.0189 0.0152 0.0096 -0.7049 0.0345 0.0168 0.0012

9.2 annualized returns

symbol AnnualizedReturn AnnualizedSharpe(Rf=0%) AnnualizedStdDev
SPY 0.0671 0.5690 0.1179
IYR 0.0152 0.1371 0.1107
QQQ 0.0837 0.5325 0.1571
TLT -0.0690 -0.9040 0.0763
GLD -0.0516 -0.7486 0.0689
IWM 0.0359 0.2328 0.1543
EEM -0.0983 -0.6160 0.1596
DBC -0.0034 -0.0300 0.1146
EFA -0.0539 -0.4514 0.1194

9.3 downside risks

symbol DownsideDeviation(0%) DownsideDeviation(MAR=10%) DownsideDeviation(Rf=0%) GainDeviation HistoricalES(95%) HistoricalVaR(95%) LossDeviation MaximumDrawdown ModifiedES(95%) ModifiedVaR(95%) SemiDeviation
SPY 0.0230 0.0269 0.0230 0.0178 -0.0691 -0.0495 0.0279 0.0706 -0.0673 -0.0543 0.0257
IYR 0.0228 0.0273 0.0228 0.0146 -0.0666 -0.0448 0.0235 0.0959 -0.0637 -0.0542 0.0237
QQQ 0.0277 0.0318 0.0277 0.0303 -0.0860 -0.0589 0.0324 0.1133 -0.0861 -0.0666 0.0315
TLT 0.0185 0.0244 0.0185 0.0085 -0.0326 -0.0313 0.0084 0.0868 -0.0430 -0.0403 0.0145
GLD 0.0152 0.0213 0.0152 0.0107 -0.0361 -0.0279 0.0109 0.1166 -0.0382 -0.0339 0.0124
IWM 0.0330 0.0367 0.0330 0.0217 -0.1099 -0.0670 0.0438 0.1339 -0.1043 -0.0782 0.0346
EEM 0.0353 0.0406 0.0353 0.0282 -0.0876 -0.0704 0.0263 0.2276 -0.0919 -0.0786 0.0306
DBC 0.0246 0.0295 0.0246 0.0117 -0.0562 -0.0547 0.0171 0.1098 -0.0613 -0.0576 0.0247
EFA 0.0279 0.0326 0.0279 0.0154 -0.0813 -0.0615 0.0275 0.1373 -0.0803 -0.0647 0.0258

9.4 downside risk ratios

symbol Annualiseddownsiderisk Downsidepotential Monthlydownsiderisk Omega Omega-sharperatio Sortinoratio Upsidepotential Upsidepotentialratio
SPY 0.0795 0.0103 0.0230 1.5770 0.5770 0.2600 0.0163 0.5693
IYR 0.0788 0.0116 0.0228 1.1497 0.1497 0.0761 0.0133 0.6734
QQQ 0.0959 0.0131 0.0277 1.5875 0.5875 0.2769 0.0207 0.7542
TLT 0.0640 0.0129 0.0185 0.5567 -0.4433 -0.3093 0.0072 0.6175
GLD 0.0526 0.0104 0.0152 0.5959 -0.4041 -0.2778 0.0062 1.1078
IWM 0.1142 0.0138 0.0330 1.2820 0.2820 0.1181 0.0177 0.5410
EEM 0.1221 0.0223 0.0353 0.6591 -0.3409 -0.2157 0.0147 0.6629
DBC 0.0852 0.0141 0.0246 1.0159 0.0159 0.0091 0.0143 0.5856
EFA 0.0967 0.0150 0.0279 0.7306 -0.2694 -0.1449 0.0110 0.4958

9.5 value at risk

symbol VaR
SPY -0.0543204
IYR -0.0541734
QQQ -0.0666155
TLT -0.0403022
GLD -0.0338613
IWM -0.0782136
EEM -0.0785758
DBC -0.0575836
EFA -0.0646908

9.6 sharpe ratio

symbol ESSharpe(Rf=0%,p=95%) StdDevSharpe(Rf=0%,p=95%) VaRSharpe(Rf=0%,p=95%)
SPY 0.0887032 0.1753105 0.1098641
IYR 0.0271849 0.0542334 0.0319830
QQQ 0.0890736 0.1690699 0.1151180
TLT -0.1327730 -0.2593925 -0.1417623
GLD -0.1103754 -0.2121981 -0.1246630
IWM 0.0373081 0.0873565 0.0497653
EEM -0.0827627 -0.1650276 -0.0967537
DBC 0.0036584 0.0067779 0.0038938
EFA -0.0503964 -0.1173963 -0.0625411

10 Appendix

10.2 Second EF Example from Stack Overflow

require(quadprog)
# min_x(-d^T x + 1/2 b^T D x) r.t A.x>=b
MV_QP <- function(nx, tarRet, Sig = NULL, long_only = FALSE) {
    if (is.null(Sig)) 
        Sig = cov(nx)
    dvec = rep(0, ncol(Sig))
    meq = 2
    Amat = rbind(rep(1, ncol(Sig)), apply(nx, 2, mean))
    bvec = c(1, tarRet)
    if (long_only) {
        meq = 1
        Amat = Amat[-1, ]
        Amat = rbind(Amat, diag(1, ncol(Sig)), rep(1, ncol(Sig)), rep(-1, ncol(Sig)))
        bvec = bvec[-1]
        bvec = c(bvec, rep(0, ncol(Sig)), 0.98, -1.02)
    }
    sol <- solve.QP(Dmat = Sig, dvec, t(Amat), bvec, meq = meq)$solution
}

steps = 50
x = stock_returns
µ.b <- apply(X = x, 2, FUN = mean)
long_only = TRUE
range.bl <- seq(from = min(µ.b), to = max(µ.b) * ifelse(long_only, 1, 1.6), 
    length.out = steps)
risk.bl <- t(sapply(range.bl, function(targetReturn) {
    w <- MV_QP(x, targetReturn, long_only = long_only)
    c(sd(x %*% w), w)
}))

weigthsl = round(risk.bl[, -1], 4)
colnames(weigthsl) = colnames(x)
weigthsl %>% head()
    DBC IWM    IYR QQQ    SPY

[1,] 0.9770 0 0.0030 0 0.0000 [2,] 0.8763 0 0.1037 0 0.0000 [3,] 0.7756 0 0.2044 0 0.0000 [4,] 0.6749 0 0.3051 0 0.0000 [5,] 0.5741 0 0.4059 0 0.0000 [6,] 0.4908 0 0.4831 0 0.0062

risk.bl = risk.bl[, 1]
rets.bl = weigthsl %*% µ.b
fan = 12
plot(x = risk.bl * fan^0.5, y = rets.bl * fan, col = 2, pch = 21, xlab = "Annualized Risk ", 
    ylab = "Annualized Return", main = "long only EF with solve.QP")

10.3 Example - Multi Layer Optimization

# demo(multi_layer_optimization)